
## ----------------------------------------------------------------- ##
## REPLICATION CODE FOR:                                             ##
## DISPROPORTIONALITY IN MEDIA REPRESENTATION OF CAMPAIGN NEGATIVITY ##
## By: Dominic Nyhuis, Hyunjin (Jin) Song, & Hajo G. Boomgaarden     ##
## ----------------------------------------------------------------- ##
## Replication Code by: Dominic Nyhuis & Hyunjin (Jin) Song          ##
## ----------------------------------------------------------------- ##

## Clear workspace, disable scientific notation
rm(list = ls())
options(scipen = 999)

## Install required packages (if not already installed)
list.of.packages <- c("stringr", "texreg", "ergm", "reshape", "xtable",
                      "data.table", "ggplot2", "ergm.count", "haven")
new.packages <- list.of.packages[!(list.of.packages %in%
                                     installed.packages()[,"Package"])]
if(length(new.packages)) install.packages(new.packages, dependencies = T)

## Automatically setting working directories
setwd(dirname(rstudioapi::getActiveDocumentContext()$path))

## Load packages
library(stringr)
library(haven)
library(parallel)
library(texreg)
library(ergm)
library(ergm.count)
library(data.table)
library(ggplot2)

## This script expects following data files:

## AUTNES 2013 media data:
# Eberl, Jakob-Moritz; Vonbun, Ramona; Haselmayer, Martin; Jacobi, Carina; Kleinen-von Königslöw, Katharina; Schönbach, Klaus; Boomgaarden, Hajo G. (2016): AUTNES Manual Content Analysis of the Media Coverage 2013. GESIS Data Archive, Cologne. ZA5864 Data file Version 1.0.0, doi:10.4232/1.12565.
# 
## AUTNES 2013 Press-release data:
# Müller, Wolfgang C.; Bodlos, Anita; Dolezal, Martin; Eder, Nikolaus; Ennser-Jedenastik, Laurenz; Kaltenegger, Matthias; Meyer, Thomas M.; Praprotnik, Katrin; Winkler, Anna Katharina (2017): AUTNES Content Analysis of Party Press Releases (OTS) 2013. GESIS Data Archive, Cologne. ZA5861 Data file Version 1.0.0, doi:10.4232/1.12688
# 
## AUTNES 2013 Party manifesto data:
# Müller, Wolfgang C.; Bodlos, Anita; Dolezal, Martin; Eder, Nikolaus; Ennser-Jedenastik, Laurenz; Kaltenegger, Matthias; Meyer, Thomas M.; Praprotnik, Katrin; Winkler, Anna Katharina (2017): AUTNES Content Analysis of Party Manifestos 2013. GESIS Data Archive, Cologne. ZA6877 Data file Version 1.1.0, doi:10.4232/1.12752
# 
## AUTNES Pre- and Post Panel Study 2013
# Kritzinger, Sylvia; Zeglovits, Eva; Aichholzer, Julian; Glantschnigg, Christian; Glinitzer, Konstantin; Johann, David; Thomas, Kathrin; Wagner, Markus (2017): AUTNES Pre- and Post Panel Study 2013. GESIS Data Archive, Cologne. ZA5859 Data file Version 2.0.1, doi:10.4232/1.12724

# DOWNLOAD AND PLACE THE STATA FILES (VERSON 14 OR HIGHER) 
# IN THE CURRENT WORKING DIRECTORY

## Set global variables
selection <- c("Budget", "Economy", "Environment", "Immigration", "Welfare")
file.names1 <- paste0("MCMC_diagnostics_valued_", selection, "%03d.pdf")
file.names2 <- paste0("GOF_valued_", selection, ".pdf")

## Data processing and load helper functions
source("Data Processing and helper functions.R")

#-----------------------
# Set up ERGM
#-----------------------

## loop over each selection, estimate model, save as lists
set.seed(12345)
fitted.models.Table1 <- lapply(1:length(selection), function(i) {
  
  media_short <- media[media$topic == selection[i],]
  ots_short <- ots[ots$topic == selection[i],]
  
  library(ergm.count)
    ## Create network
    media_mat <- as.matrix(table(media_short$subject, media_short$object))
    ## Discard loops
    diag(media_mat) <- 0
  
    ## Create statnet network and set count edge
    g <- as.network(
      media_mat, 
      directed = T, 
      ignore.eval = FALSE, 
      names.eval = "edge_count"
    )
  
    ## Set up press releases network
    ots_mat <- as.matrix(table(ots_short$subject, ots_short$object))
    diag(ots_mat) <- 0
    colnames(ots_mat) <- rownames(ots_mat) <- rownames(media_mat)
    class(ots_mat) <- "matrix"
      
    coalitions <- ots_mat
    coalitions[coalitions > 0] <- 0
    coalitions[1,2] <- 1
    coalitions[2,1] <- 1

    ## incomvency indicators
    g %v% "chancellor" <- c(1,0,0,0,0,0)
    g %v% "government" <- c(1,1,0,0,0,0)
      
    ## Electoral result
    g %v% "party_size" <- c(26.8, 24.0, 20.5,
                            12.4, 3.5, 5.7)/100

    ## manifesto-based issue importance
    g %v% "emphasis"  <- issue_emphasis[,selection[i]]
    ## competence based issue ownership
    g %v% "competence_prop" <- competence_prop[,selection[i]]
    ## competence based issue ownership
    g %v% "competence_bin" <- competence_bin[,selection[i]]
    ## party ideological placement
    g %v% "lr" <- lr

#-----------------------
# Model Table 1
#-----------------------
    RNGkind("L'Ecuyer-CMRG")
    set.seed(12345)
    
  fit <- ergm(
    g ~ sum + # baseline
      nodeosqrtcovar + edgecov(ots_mat, attrname = "pressreleases", form = "sum") +
      nodeifactor("competence_bin", form = "sum") + ## perceived issue ownership (binary)
      nodeofactor("competence_bin", form = "sum") + ## perceived issue ownership (binary)
      nodeifactor("government", form = "sum") + 
      nodeocov("party_size", form = "sum") + absdiff("lr")
    ,
    estimate = "MLE", 
    response = "edge_count",
    reference = ~ Poisson, 
    verbose = FALSE, 
    control = control.ergm(
      MCMLE.trustregion = 200,
      MCMC.burnin= 1024 * 2 * 16,
      MCMC.interval= 1024 * 2,
      MCMLE.conv.min.pval = 0.05
    )
  )

  return.list <- list(fit = fit, logLik = fit$mle.lik, ots_mat = ots_mat)
  return(return.list)
})

## Print MCMC diagnostic plots (Figure S1 to S5 in online SI)
for (i in 1:5) {
  file.name <- file.names1[i]
  ## Model goodness of fit
  pdf(file = file.name, width = 5.8, height = 8.3, onefile = F)
  mcmc.diagnostics(fitted.models.Table1[[i]]$fit, 
                   vars.per.page = 4, verbose = F)
  dev.off()
}

## manually plot GOF against simulated networks, 
## using additional statistics that are not explicitly modeled in the fit object
## (Figure S6 to S10 in online SI)
require(ggplot2)
for (i in 1:5) {
  fit <- fitted.models.Table1[[i]]$fit
  g <- fitted.models.Table1[[i]]$fit$network
  ots_mat <- fitted.models.Table1[[i]]$ots_mat
  
  gof.fit <- gof.valued.ergm(fit)
    pp <- ggplot() + stat_boxplot(data = gof.fit$estsim, aes(x = terms, y = y), 
                                  geom = "errorbar", width = 0.5, color = "gray50") + 
      geom_boxplot(data = gof.fit$estsim, aes(x = terms, y = y), 
                   outlier.shape = NA, color = "gray50") +
      theme_bw() + ylab("Simulated Quantiles (centered)") + 
      xlab("") + geom_point(data = gof.fit$estobs, aes(x = terms, y = y), size = 2) + 
      geom_line(data = gof.fit$estobs, aes(x = terms, y = y, group = 1), size = 1)
    
    file.name <- file.names2[i]
    ggsave(file = file.name)

}

## Table 1 in the manuscript
texreg::screenreg(
  sapply(fitted.models.Table1, "[", 1),
  custom.model.names = selection,
  custom.coef.names = c("Sum", "Individual heterogeneity", 
                        "Objective attacks", "Issue owner (binary, Target)",
                        "Issue owner (binary, Attacker)", "Incumbent (Target)", 
                        "Party size", "Ideological distance"),
  dcolumn = T, booktabs = T
)


#-----------------------
# Online SI: 
#-----------------------

## Descriptive statistics by issue area (Table S1)
descriptive <- rbind(
  ## Ties
  sapply(1:length(selection), function(i) {
    media_short <- media[media$topic == selection[i],]
    media_mat <- as.matrix(table(media_short$subject, media_short$object))
    diag(media_mat) <- 0
    sum(media_mat)}),
  ## Unique Ties
  sapply(1:length(selection), function(i) {
    media_short <- media[media$topic == selection[i],]
    media_mat <- as.matrix(table(media_short$subject, media_short$object))
    diag(media_mat) <- 0
    media_mat <- sna::event2dichot(media_mat,method = "absolute", thresh = 0)
    sum(media_mat)}),
  ## Average no. of attacks, sd, max and min
  sapply(1:length(selection), function(i) {
    media_short <- media[media$topic == selection[i],]
    media_mat <- as.matrix(table(media_short$subject, media_short$object))
    diag(media_mat) <- 0
    return(c("mean" = mean(media_mat), "sd" = sd(as.vector(media_mat)), 
             "max" = max(as.vector(media_mat)), "min" = min(as.vector(media_mat))))
  })
)

colnames(descriptive) <- selection
rownames(descriptive) <- c("Ties","Ties (unique)", "Average", "SD", "Max", "Min")
descriptive ## Appendix Table S1 

## Newspapers
## 3101 (Der Standard)  
## 3102 (Die Presse + Die Presse am Sonntag) 
## 3103 (Salzburger Nachrichten)  
## 3201 (Kronen Zeitung)  
## 3202 (Österreich)  
## 3203 (Heute)  
## 3301 (Kurier)  
## 3401 (Kleine Zeitung)
np.list <- c(3101, 3102, 3103, 3201, 3202, 3203, 3301, 3401)
descriptive_media <- rbind(
  ## Ties
  sapply(1:length(np.list), function(i) {
    media_short <- media[media$media == np.list[i],]
    media_mat <- as.matrix(table(media_short$subject, media_short$object))
    diag(media_mat) <- 0
    sum(media_mat)}),
  ## Unique Ties
  sapply(1:length(np.list), function(i) {
    media_short <- media[media$media == np.list[i],]
    media_mat <- as.matrix(table(media_short$subject, media_short$object))
    diag(media_mat) <- 0
    media_mat <- sna::event2dichot(media_mat,method = "absolute", thresh = 0)
    sum(media_mat)}),
  ## Average no. of attacks, sd, max and min
  sapply(1:length(np.list), function(i) {
    media_short <- media[media$media == np.list[i],]
    media_mat <- as.matrix(table(media_short$subject, media_short$object))
    diag(media_mat) <- 0
    return(c("mean" = mean(media_mat), "sd" = sd(as.vector(media_mat)), 
             "max" = max(as.vector(media_mat)), "min" = min(as.vector(media_mat))))
  })
)

rownames(descriptive_media) <- c("Ties","Ties (unique)", "Average", "SD", "Max", "Min")  
colnames(descriptive_media) <- c("DS", "DP", "SN", ## Quality
                                 "Kronen", "Ö", "Heute", ## Tabloid
                                 "Kurier", "Kleine") ## mid-range paper

## Descriptive statistics by newspapers (Table S2)
descriptive_media

#-----------------------
# Table S3: proportion measure of ownership
#-----------------------

## loop over each selection, estimate model, save as lists
fitted.models.TableS3 <- lapply(1:length(selection), function(i) {
  
  assign("g", fitted.models.Table1[[i]]$fit$network)
  assign("ots_mat", fitted.models.Table1[[i]]$ots_mat)
  RNGkind("L'Ecuyer-CMRG")
  set.seed(578923473)
fit <- ergm(
  g ~ 
    sum + # baseline
    nodeosqrtcovar +
    edgecov(ots_mat, attrname = "pressreleases", form = "sum") +
    nodeocov("competence_prop", form = "sum") + 
    nodeicov("competence_prop", form = "sum") + 
    nodeifactor("government", form = "sum") +
    nodeocov("party_size", form = "sum") +
    absdiff("lr", form = "sum")
  ,
  response = "edge_count",
  reference = ~ Poisson, 
  verbose = FALSE, 
  #estimate = "MLE", 
  control = control.ergm(
    MCMLE.trustregion = 200,
    MCMC.burnin= 1024 * 2 * 16,
    MCMC.interval= 1024 * 2,
    CD.conv.min.pval = 0.05
  )
)

  return.list <- list(fit = fit, logLik = fit$mle.lik, ots_mat = ots_mat)
  return(return.list)
})

## Table S3, proportion measure
texreg::screenreg(
  sapply(fitted.models.TableS3, "[", 1),
  custom.model.names = selection,
  custom.coef.names = c("Sum", "Individual heterogeneity", 
                        "Objective attacks", "Issue owner (Prop., Target)",
                        "Issue owner (Prop., Attacker)", "Incumbent (Target)", 
                        "Party size", "Ideological distance")
)

#-----------------------
# Model Table S4: baseline control over/under dispersion
#-----------------------

set.seed(578923473)
fitted.models.TableS4 <- lapply(1:length(selection), function(i) {
  
  assign("g", fitted.models.Table1[[i]]$fit$network)
  assign("ots_mat", fitted.models.Table1[[i]]$ots_mat)
  RNGkind("L'Ecuyer-CMRG")
  set.seed(578923473)
  fit <- ergm(
    g ~ 
      sum(pow = 1/2) + # baseline control over/under dispersion
      nodeosqrtcovar +
      edgecov(ots_mat, attrname = "pressreleases", form = "sum") +
      nodeifactor("competence_bin", form = "sum") + ## perceived issue ownership (binary)
      nodeofactor("competence_bin", form = "sum") + ## perceived issue ownership (binary)
      nodeifactor("government", form = "sum") +
      nodeocov("party_size", form = "sum") +
      absdiff("lr")
    ,
    response = "edge_count", reference = ~ Poisson, verbose = FALSE, estimate = "MLE", 
    control = control.ergm(
      MCMLE.trustregion = 200,
      MCMC.burnin= 1024 * 2 * 16,
      MCMC.interval= 1024 * 2,
      MCMLE.conv.min.pval = 0.05
    )
  )
  
  return.list <- list(fit = fit, logLik = fit$mle.lik, ots_mat = ots_mat)
  return(return.list)
})

texreg::screenreg(
  sapply(fitted.models.TableS4, "[", 1),
  custom.model.names = selection,
  custom.coef.names = c("Sum (fractional moment)", "Individual heterogeneity", 
                        "Objective attacks", "Issue owner (Target)",
                        "Issue owner (Attacker)", "Incumbent (Target)", 
                        "Party size", "Ideological distance"),
  dcolumn = T, booktabs = T
)




## ------------------------------------##
## Table A5: Jackknifing Media Outlets
## ------------------------------------##

## Newspapers
## 3101 (Der Standard)  
## 3102 (Die Presse + Die Presse am Sonntag) 
## 3103 (Salzburger Nachrichten)  
## 3201 (Kronen Zeitung)  
## 3202 (Österreich)  
## 3203 (Heute)  
## 3301 (Kurier)  
## 3401 (Kleine Zeitung)
np.list <- list(3101, 3102, 3103, 3201, 3202, 3203, 3301, 3401)

require(parallel)
jackknifing.ergm <- mclapply(np.list, function(x) {
  
  ## loop over each selection, estimate model, save as lists
  fitted.models.Table1.jackknifing <- lapply(1:length(selection), function(i) {
    
    media_short <- media[(media$topic == selection[i]) & (media$media != x),]

    library(ergm.count)
    ## Create network
    media_mat <- as.matrix(table(media_short$subject, media_short$object))
    diag(media_mat) <- 0
    
    ## Create statnet network and set count edge
    g <- as.network(
      media_mat, 
      directed = T, 
      ignore.eval = FALSE, 
      names.eval = "edge_count"
    )
    
    ## incomvency indicators
    g %v% "chancellor" <- c(1,0,0,0,0,0)
    g %v% "government" <- c(1,1,0,0,0,0)
    
    ## Electoral result
    g %v% "party_size" <- c(26.8, 24.0, 20.5,
                            12.4, 3.5, 5.7)/100
    
    ## manifesto-based issue importance
    g %v% "emphasis"  <- issue_emphasis[,selection[i]]
    ## competence based issue ownership
    g %v% "competence_prop" <- competence_prop[,selection[i]]
    ## competence based issue ownership
    g %v% "competence_bin" <- competence_bin[,selection[i]]
    ## party ideological placement
    g %v% "lr" <- lr
    
    ## press release
    assign("ots_mat", fitted.models.Table1[[i]]$ots_mat)
    
    RNGkind("L'Ecuyer-CMRG")
    set.seed(578923473)
    
    fit <- ergm(
      g ~ 
        sum + # baseline
        nodeosqrtcovar +
        edgecov(ots_mat, attrname = "pressreleases") +
        nodeifactor("competence_bin") + ## perceived issue ownership (binary)
        nodeofactor("competence_bin") + ## perceived issue ownership (binary)
        nodeifactor("government") +
        nodeocov("party_size") +
        absdiff("lr")
      ,
      response = "edge_count",
      reference = ~ Poisson, 
      verbose = FALSE, 
      control = control.ergm(
        MCMLE.trustregion = 200,
        MCMC.burnin= 1024 * 2 * 16,
        MCMC.interval= 1024 * 2,
        MCMLE.conv.min.pval = 0.05
      )
    )
    
    return.list <- list(fit = fit, logLik = fit$mle.lik, ots_mat = ots_mat)
    return(return.list)
  })
  
  ergm.coef <- lapply(fitted.models.Table1.jackknifing, function(y) {
    ergm.output <- summary(y$fit)$coefficient[, c(1,2,5)]
    ergm.output
  })
  ergm.coef <- do.call(rbind, ergm.coef)  
  colnames(ergm.coef) <- c("est", "se", "p")
  ergm.coef$topic <- rep(selection, 8)
  
  ergm.gof <- lapply(fitted.models.Table1.jackknifing, function(y) {
    aic <- summary(y$fit)$aic
    bic <- summary(y$fit)$bic
    lik <- y$fit$mle.lik[1]
  ergm.gof <- c(aic = aic, bic = bic, lik = lik)
  ergm.gof})
  
  return.list <- list(ergm.coef = ergm.coef, ergm.gof = ergm.gof)  
  return(return.list)
}, 
  mc.cores = parallel::detectCores()) 

## process output
jackknifing.ergm.gof <- do.call(rbind, 
                                unlist(sapply(jackknifing.ergm, "[", 2), 
                                       recursive = F)) %>% as.data.frame(.)
jackknifing.ergm.gof$topic <- rep(selection, 8)
setDT(jackknifing.ergm.gof)

jackknifing.ergm <- do.call(rbind, sapply(jackknifing.ergm, "[", 1)) %>% as.data.frame(.)
jackknifing.ergm$vars <- rep(names(coef(fitted.models.Table1[[1]]$fit)), times = 8)
jackknifing.ergm$omitted.outlet <- rep(unlist(np.list), each = 8)
setDT(jackknifing.ergm)

require(texreg)
fitted.models.TableS5 <- make.texreg.jackknifing.models()
  
## comparison of results
screenreg(fitted.models.TableS5, ## jacknifing outlet models
          custom.model.names = selection,
          custom.coef.names = c("Sum", "Individual heterogeneity", 
                                "Objective attacks", "Issue owner (Binary, Target)",
                                "Issue owner (Binary, Attacker)", 
                                "Incumbent (Target)", 
                                "Party size", "Ideological distance")) 

## --------------------------------------------------- ##
## Table S6: Models excluding ideological distance     
## --------------------------------------------------- ##

fitted.models.TableS6 <- lapply(1:length(selection), function(i) {
  
  assign("g", fitted.models.Table1[[i]]$fit$network)
  assign("ots_mat", fitted.models.Table1[[i]]$ots_mat)
  
  RNGkind("L'Ecuyer-CMRG")
  set.seed(578923473)
  fit <- ergm(
    g ~ 
      sum + # baseline
      nodeosqrtcovar +
      edgecov(ots_mat, attrname = "pressreleases") +
      nodeocov("competence_bin") + ## perceived issue ownership (proportion)
      nodeicov("competence_bin") + ## perceived issue ownership (proportion)
      nodeifactor("government") +
      nodeocov("party_size")
    ,
    response = "edge_count",
    reference = ~ Poisson, 
    verbose = FALSE, estimate = "MLE", 
    control = control.ergm(
      MCMLE.trustregion = 200,
      MCMC.burnin= 1024 * 2 * 16,
      MCMC.interval= 1024 * 2,
      MCMLE.conv.min.pval = 0.05
    )
  )
  
  return.list <- list(fit = fit, logLik = fit$mle.lik, ots_mat = ots_mat)
  return(return.list)
})

screenreg(sapply(fitted.models.TableS6, "[", 1),
          custom.model.names = selection,
          custom.coef.names = c("Sum", "Individual heterogeneity", 
                                "Objective attacks", "Issue owner (Binary, Target)",
                                "Issue owner (Binary, Attacker)", 
                                "Incumbent (Target)", 
                                "Party size"))


## Visulaization of networks (figure S11)
library(igraph)
pdf(file = "network.plots.pdf", paper = "a4")
par(cex = 0.7, mai = c(0.1,0.1,0.1,0.1))
lapply(1:length(selection), function(i) {
  media_short <- media[media$topic == selection[i],]
  media_mat <- as.matrix(table(media_short$subject, media_short$object))
  g <- igraph::graph_from_adjacency_matrix(media_mat, weighted = T)
  g <- simplify(g, edge.attr.comb = "sum")
  V(g)$size <- c(26.8, 24.0, 20.5, 12.4, 3.5, 5.7)
  V(g)$color <- c("#ce000c", "#000000", "#0E428E", "#6BA325", "#F29400", "#F8D323")
  V(g)$label.color <- "black"
  plot(g, edge.width = 1 + edge_attr(g)$weight/10, layout = layout_in_circle,
          edge.arrow.size = 0.7, vertex.label.dist = 3, edge.curved = .15,
       main = paste0("Issue area: ", selection[i]))
})
  media_short <- media[media$topic %in% selection,]
  media_mat <- as.matrix(table(media_short$subject, media_short$object))
  g <- igraph::graph_from_adjacency_matrix(media_mat, weighted = T)
  g <- simplify(g, edge.attr.comb = "sum")
  V(g)$size <- c(26.8, 24.0, 20.5, 12.4, 3.5, 5.7)
  V(g)$color <- c("#ce000c", "#000000", "#0E428E", "#6BA325", "#F29400", "#F8D323")
  V(g)$label.color <- "black"
  plot(g, edge.width = 1 + edge_attr(g)$weight/10, layout = layout_in_circle,
       edge.arrow.size = 0.7, vertex.label.dist = 3, edge.curved = .15,
       main = "Issue area: Total")
dev.off()

  

 
 





















